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Hotspots of genetic diversity are regions of utmost importance for species survival and conservation, and 
their intimate link with the geographic location of glacial refugia has been well established. Nonetheless, the 
microevolutionary processes underlying the generation of hotspots in such regions have only recently 
become a fervent field of research. We investigated the phylogeographic and population genetic structure of 
the agile frog, Rana dalmatina, within its putative refugium in peninsular Italy. We found this region to 
harbour far more diversity, phylogeographic structure, and lineages of ancient origin than that by the rest of 
the species' range in Europe. This pattern appeared to be well explained by climate-driven 
microevolutionary processes that occurred during both glacial and interglacial epochs. Therefore, the 
inferred evolutionary history of R. dalmatina in Italy supports a view of glacial refugia as 'factories' rather 
than as repositories of genetic diversity, with significant implications for conservation strategies for 
hotspots. 

Hotspots of genetic diversity, defined as geographic areas harbouring a major portion of a species' genetic 
diversity, have been observed for most species studied to date. The discovery of hotspots and the causes of 
their formation have become vigorous areas of research 110 . Indeed, hotspots are increasingly seen as 
excellent opportunities for the study of long-debated themes in ecology and evolutionary biology, including the 
microevolutionary processes leading to the formation of latitudinal patterns of biodiversity, the role of divergence 
with or without gene flow in moulding current biodiversity, the role of hybridization in evolution, and the 
geographic mosaic of revolutionary interactions 5,6 ' 1,12 . On the other hand, these areas are also of the utmost 
importance for long-term survival and conservation of species, since they frequently harbour most species' 
resources for coping with natural and anthropogenic environmental changes 1315 . 

Two main evolutionary scenarios have been described to date for the generation of hotspots of intraspecific 
biodiversity in temperate species 16 . Under each of these scenarios, past climatic oscillations are seen as the major 
distal cause of the process and the location of hotspots is seen as intimately linked to the so-called 'glacial refugia'. 
Under the first scenario, a hotspot would testify to the long-term persistence of demographically stable, large 
populations within a glacial refugium. In contrast, under the second scenario the high genetic diversity found in 
glacial refugia would be the outcome of climate-driven cycles of allopatric differentiation within sub-refugia, 
which is possibly followed by a secondary admixture of divergent lineages. In spite of some fossil evidence 
suggesting that refugial populations could have been small, isolated, and at low density 1718 , the first scenario 
was considered as prevalent in early literature, particularly for wide-ranging species 8,19,20 . However, mounting 
data from phylogeographic studies focused on populations within glacial refugia identify the second scenario as 
relevant for many temperate species and geographic regions 1,5,7,9,21,22 . 

Identifying the geographic locations of hotspots of genetic diversity, the genetic structure of populations 
therein, and the microevolutionary processes that have led to their formation are important not only for the 
advancement of knowledge in ecology, biogeography, and evolutionary biology, but also have major practical 
implications. For instance, regional networks of protected areas and management plans for species conservation 
are often designed under the assumption of homogeneous distribution of diversity within species' populations, 
even when there is inadequate knowledge of the geography and depth of the population genetic structure 14,15 . Such 
regional networks and plans would be appropriate only in scenarios where a single evolutionary unit (ESU; 
sensu 23 ) is present, whereas they should be carefully re-evaluated in cases where multiple ESUs are found, 
deserving separate conservation efforts, or in cases where major hotspots of genetic diversity fall outside such 
networks. 
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Table 1 | Geographic location of the 40 sampled populations of Rana dalmatina, number of individuals analysed using mtDNA(n mt ) and 
microsatellite (n ms ) genetic markers, and distribution of the 49 identified mtDNA haplotypes among populations 



Sampl 






Latitude N 


Longitude E 


rirYit/ rime 

1 1 ml/ 1 1 ms 


Haplotypes (n) 


1 


Slovenia 


Komen 


45=49' 


13=45' 


3/3 


Nil (3) 


2 




Gorjansko 


45=48' 


13=42' 


1/1 


Nll(l) 


3 




Kocevje 


45=37' 


14=53' 


1/1 


Mil 1(1) 


4 


Croazia 


Roc 


45=23' 


14 03' 


1/1 


Nll(l) 


5 




Livade-Trombal 


45=21' 


13=47' 


3/1 


Nil (3) 


6 




Svetvincenant 


45=04' 


13=53' 


4/4 


Nil (3), NI8(1) 


7 




Dreznica 


45=06' 


15=07' 


4/4 


Nll(l), NI9(1), Nil 0(1 ), Nil 2(1 ) 


8 


Italia 


Punta Alberete 


44=31' 


12=13' 


6/6 


NI1(5),NI5(1) 


9 




Pianoro 


44=29' 


1 1=20' 


4/4 


NI1(3),NI6(1) 


10 




Pineta di Classe 


44° 2 r 


12=19' 


6/6 


Nil (5), NI7(1) 


1 1 




Greve in Chianti 


43=35' 


11=18' 


4/4 


Nil (4) 


12 




Torrente Farma 


43° 47' 


11=14' 


3/3 


Nil (3) 


13 




Trasimeno 


43=08' 


12=10' 


1/1 


Mil (1 ) 


14 




Urbino 


43=43' 


12 38' 


1/1 


NI4(1) 


15 




Selva del Lamone 


42=33' 


11=43' 


7/- 


Nil (4), Nil 3(3) 


16 




Canale Monterano 


42=08' 


12=06' 


3/3 


Nil (3) 


17 




Stagni delta Doganella 


41=42' 


12=43' 


8/5 


Nil (7), Nlll(l) 


18 




Decima Malafede 


41=44' 


12=25' 


6/3 


Nil (4), NI2(1), NI3(1), 


19 




Foglino 


41=27' 


12=42' 


8/8 


Nll(l), Nil 1 (4), NII2(2), NII3(1) 


20 




Circeo 


41°14' 


13=04' 


2/2 


NII2(2) 


21 




Pescolanciano 


41=40' 


14=20' 


1/1 


SI8(1) 


22 




Foresta Umbra 


41=51' 


16=01' 


8/8 


SI2(6), SI3(1), SI4(1) 


23 




Lago Fondo 


39 55' 


16=14' 


8/5 


SI5(5), SI6(2), SI7(1) 


24 




Acqualisparti 


39=59' 


15=52' 


6/5 


SI6(2), SI7(4) 


25 




Massadita 


39=51' 


16=18' 


6/5 


SI6(6) 


26 




Santa Domenica Talao 


39=49' 


15=52' 


8/6 


Sll(l), SI9(7) 


27 




Orsomarso 


39=48' 


15=54' 


1/- 


Sll(l) 


28 




Ficara 


39=39' 


16=07' 


4/4 


SI9(1), SII9(2), SII15(1) 


29 




Lago Trifoglietti 


39=33' 


16=01' 


9/9 


SII5(1), SII6(2), SII7(2), SII8(1), SII12(1), Sill 3(2) 


30 




San Lucido 


39=17' 


16=05' 


19/5 


SI10(1), SII4(1), SII9(1), Sill 3(2), SIM 4(3), SIM 8(4), SII22(7) 


31 




Cerisano 


39=15' 


1 6=09' 


8/8 


SII5(2), Sill 8(3), 51122(3) 


32 




Belmonte Calabro 


39=10' 


16=05' 


13/5 


SI9(2), Sill (6), SII20(2), 51122(3) 


33 




San Pietro in Guarano 


39=20 


16=18' 


7/5 


SII5(2), Sll 1 3(2), SII19(2), 51122(1 ) 


34 




Campo San Lorenzo 


39=20' 


16=29' 


1 1/6 


SII5(6), SII9(1), SII19(1), SII22(3) 


35 




Macchialonga 


39=22' 


16=32' 


12/4 


SII5(1 1), SII7(1) 


36 




Lago di Angitola 


39=07' 


16=06' 


6/6 


Slll(2), SII9(1), SIM 1(2), SII16(1) 


37 




Capo Vaticano 


38=39' 


16=01' 


10/7 


Sill (4), SII9(5), SII23(1) 


38 




Serra San Bruno 


38=36' 


16=20' 


18/8 


Sill (7), SII3(3), SIM 0(1), Sill 7(2), SII9(2), Sill 6(2), SII21 (1) 


39 




Zomaro 


38=40' 


15=59' 


3/3 


Sill (3) 


40 




Gambarie 


38=10' 


15=50' 


10/9 


Sill (9), SII2(1) 



The agile frog, Rana dalmatina Fitzinger in Bonaparte, 1839, is a 
brown frog that is widely distributed in western, central, and south- 
ern Europe, where it inhabits a variety of landscapes, mostly within 
mixed deciduous forests from sea level to 2000 m a.s.l. R. dalmatina 
spends most of its life on land, breeding in a variety of natural and 
anthropogenic water bodies, including flooded meadows, ponds, 
lakes, temporary pools, and slow-running waters 24 . Although this 
species is not considered to be under imminent conservation risk, 
population declines and disappearances have recently been 
observed 25 . In a recent survey of the genetic variation throughout 
the species' range, R. dalmatina populations have appeared substan- 
tially homogeneous, with the single exception of a strongly differen- 
tiated population in southern Italy, which was considered as evidence 
for a glacial refugium of this species being located within this area 26 . 

In this study, we further investigate the population genetic struc- 
ture and evolutionary history of -R. dalmatina, focusing on its putat- 
ive refugial range in peninsular Italy. This geographic area has 
recently emerged as a hotspot of genetic diversity and as a site of 
multiple refugia for many temperate animal species, including 
amphibians 6-21,27-31 . Here we adopted a dense sampling strategy 
within this area, with the following aims: i) to assess the phylogeo- 
graphic and population genetic structure, and the historical demo- 
graphic trends of-R. dalmatina within its prospected glacial refugium, 
ii) to investigate the microevolutionary processes that have moulded 



its genetic diversity within this area; and iii) to evaluate the main 
implications of our results for the conservation of this species. 

Results 

Mitochondrial DNA. The final dataset included 741 bp of the 
mitochondrial cytochrome oxidase subunit I gene (COI; 46 
variable positions, 37 parsimony informative) and 653 bp for the 
mitochondrial cytochrome B gene (CYTB; 55 variable positions, 42 
parsimony informative), for all the 244 individuals analysed. Forty- 
nine composite haplotypes were found (see Table 1), showing a 
sequence divergence (p-distance) ranging between 0.1% and 4.5%. 

The statistical parsimony analysis did not connect all haplotypes 
into a single phylogenetic network (Figure 1A). Instead, 2 main 
groups of haplotypes were observed (haplogroups N and S) having 
an ML-corrected net sequence divergence of 4.0% (0.5% S.E.), each 
showing a discrete geographic distribution. Haplogroup N was found 
in north-central and northern Italy, Croatia, and Slovenia, (samples 
1-19), whereas haplogroup S was distributed in southern Italy (sam- 
ples 20-40). Both of these haplogroups were further divided into 2 
sub-haplogroups, the distribution of which is presented in Figure 1C. 
While haplogroups N and S were never found in syntopy, the respect- 
ive sub-haplogroups (NI-NII and SI-SII) showed overlapping distri- 
butions mostly limited to geographically intermediate populations. 
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Figure 1 | (A) Phylogenetic networks of the 49 mtDNA haplotypes found among Rana dalmatina populations in peninsular Italy, based on the 
statistical parsimony procedure implemented in TCS. Circle sizes are proportional to haplotype frequency (see inset, lower left); missing 
intermediate haplotypes are shown as open dots. (B) Haplotype genealogy yielded using HaplotypeViewer, based on the ML phylogenetic tree of the 49 
haplotypes found. (C) Geographic locations of the 40 sampled populations of R. dalmatina and frequency distribution of the main haplogroups 
among populations, shown as pie diagrams. Populations are numbered as in Table 1. The map was drawn using the software Canvas 1 1 (ACD Systems of 
America, Inc.). 



The haplotype genealogy yielded by HaplotypeViewer showed the 
same groupings as the statistical parsimony networks, with minor 
topological differences that did not affect the general pattern 
(Figure IB). Within this genealogy, haplogroups N and S were sepa- 
rated by 49 mutational steps. 

Estimates of genetic diversity for the 4 haplogroups are shown in 
Table 2. Particularly low values of h and n were observed for the 
northernmost haplogroup (Nl), whereas high values of both indices 
were found for the southernmost haplogroup (SII), and to a lesser 
extent, for SI. 

The AMOVA analysis was carried out by separating populations 
according to the geographic and frequency distribution of the 4 main 
haplogroups ([1-18], [19-20], [21-27], [28-40]). This analysis 
indicated that the amount of genetic variance explained by the 



among-grouplevelofvariationwas91.6% (FCT = 0.92), the amount 
explained by the variation among populations within groups was 
2.7% (FSC = 0.32), while the within-population level explained 
5.7% (FST = 0.94). All variance components and fixation indices 
were highly significant (P < 0.001). 

The analysis of the divergence time dating with BEAST yielded 
effective sample size (ESS) values that largely exceeded 200 for 
all parameters of interest and full convergence for parameter esti- 
mates between runs. The divergence time between haplogroups N 
and S was estimated at 2.109 My (95% HPD: 1.433-2.929 My). 
On the other hand, the divergence between haplogroups NI-NII 
and SI— SII was estimated to have occurred at 0.381 My (95% 
HPD: 0.172-0.594 My) and 0.346 My (95% HPD: 0.197- 
0.519 My), respectively. 



Table 2 | Estimates of haplotype (h) and nucleotide [n) diversity 60 and neutrality test statistics F/ 4 , D 76 and R2 5 , for the 4 main mtDNA 
haplogroups observed among the 244 Rana dalmatina specimens analysed 



Haplogroup 


n 


h 


n 


Fs 


D 


R 2 


Nl 


66 


0.380 (0.077) 


0.00032 (0.00032) 


-16.488** 


-2.319** 


0.031* 


Nil 


10 


0.644 (0.101) 


0.00094 (0.00072) 


1.021 


0.850 


0.21 1 


SI 


42 


0.854 (0.026) 


0.00150 (0.00095) 


-2.103 


-0.540 


0.093 


SII 


126 


0.879 (0.016) 


0.00246 (0.00140) 


-6.385 


-0.839 


0.065 


n, number of ind 
*P < 0.05; 


viduals carrying haplotypes from a 


specific haplogroup; 











**P < 0.01 . The F s statistic was considered statistically significant at a < 0.05 under nominal significance values P < 0.02 7 ' 
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Figure 2 | Bayesian skyline plots showing the historical demographic trends of 3 of the main mtDNA haplogroups identified within Rana dalmatina. 

The x-axis is scaled in years (bp). The y-axis is scaled in units of effective population size per generation length. 



The historical demographic trends for all the main haplogroups 
were investigated by means of a bayesian skyline plot (BSP) analysis, 
except for Nil because of its low sample size (n = 10). The BSPs 
estimated for the remaining haplogroups are shown in Figure 2. The 
time to the most recent common ancestors (TMRCAs) was estimated 
as 23 000 years (95% HPD: 8300-45 600), 84 500 years (95% HPD: 
27 700-158 800), and 155 000 years (95% HPD: 63 200-266400) for 
haplogroups NI, SI, and SII, respectively. For haplogroup NI, a 
marked demographic expansion was inferred that began shortly after 
the TMRCA and lasted until recent times. In contrast, the BSPs for 
haplogroups SI and SII showed no marked demographic changes; the 
median estimates of effective population size for these groups never 
fell outside the 95% HPDs estimated along the entire time axis. 

The estimated values of the neutrality test statistics are given in 
Table 2. In agreement with the BSP analysis, large negative values of 



F s and D and a small positive value of R 2 (all suggesting recent and 
marked population growth) were only observed for haplogroup NI 
(all P < 0.05). On the other hand, none of the estimated values of the 
statistics for the haplogroups Nil, SI, and SII was significant. 

Microsatellites. A total of 170 individuals from 38 locations were 
genotyped at 6 microsatellite loci (see Table 1). Among these 
individuals, none was scored at less than 5 loci, and the overall 
proportion of missing data in the final dataset was 0.2%. 

Evidence of Hardy- Weinberg or linkage disequilibria was not 
found across samples or loci. For all populations, the observed het- 
erozygosity ranged from 0.69 (HI locus) to 0.89 (G12 locus), the 
expected heterozygosity ranged from 0.84 (HI locus) to 0.96 (E8 
locus), and the number of alleles observed at each locus was as 
follows: F5 = 19, G12 = 28, Gil = 16, E8 = 38, HI = 16, and B5 
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Figure 3 | Genetic structure of Rana dalmatina populations in Italy, which were estimated using the Bayesian clustering method implemented in TESS. 

(A) Mean values of the DIC statistics (averaged over 100 runs) estimated for models with the number of genetic clusters (K) ranging from 2 to 7; 

(B) Posterior predictive map of the admixture proportions generated from the spatial interpolation procedure implemented in TESS; (C) Admixture 
proportions of each sampled individual for the 4 recovered genetic clusters. The map was drawn using the software TESS Ad-Mixer. 



= 22. Consequently, the total number of alleles found within the 
study area (n = 139) was 12.9% higher than that previously reported 
for the same loci throughout the species range 26 . 

The analysis of population genetic structure by using TESS indi- 
cated K = 4 as the best grouping option for our data, since only a 
minor decrease in DIC values was observed at higher K values 
(Figure 3A). As previously observed for the mtDNA data, 2 groups 
were found in south-central and southern Italy (samples 21-40; see 
Figure 3B). However, unlike the mtDNA data, the 2 population 
groups were largely admixed in most of the areas (Figures 3B and 
3C), with some exceptions in marginal samples. Two groups were 
observed in the north-central and northern portion of the study area. 
However, the geographic distribution of the 2 groups did not con- 
form to the geographic pattern observed in the mtDNA data. Indeed, 
the area of admixture between the 2 groups was not located in central 
Italy (samples 17-20), but on the north-eastern side of the Italian 
Peninsula (samples 8-10, 14). Finally, some instances of unidir- 



ectional gene exchange from south to north were observed among 
samples located in central Italy (samples 17-20), which otherwise 
belonged to the north-central group (Figure 3C). Finally, the mean 
expected heterozygosity estimated for the 4 groups was 0.90 (0.01 
SE), 0.92 (0.01 SE), 0.88 (0.03 SE), 0.88 (0.02 SE) for the southern, 
south-central, north-central and northern groups respectively. 

Discussion 

Our analysis of genetic variation in the agile frog in peninsular Italy 
allowed us to identify this area as not only a site of the most divergent 
lineage but also as the major hotspot of genetic diversity within this 
species. Indeed, populations from this area harbour far more genetic 
variation than that found in the central, eastern, and southern 
European populations combined 26 . Further, most of the genetic 
variation was accounted for by differences between 4 genetically 
divergent and geographically well-defined lineages, suggesting a 
long-lasting and multi-episode evolutionary history for _R. dalmatina 
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within this area. In turn, the occurrence of 4 lineages with relatively 
ancient divergence clearly indicates that peninsular Italy did not act 
as a single refugium for the species, but rather as a mosaic of refugia, 
at least during the last Pleistocene glacial-interglacial cycles. 

The most ancient event inferred in the history of R. dalmatina 
populations was the divergence between the N and S mtDNA hap- 
logroups, which has been roughly dated at 2.1 million years (MY). 
The precise dating of this divergence was prevented by the large 
95% CI for the time estimate (1.4-2.9 My), our use of an exter- 
nally calibrated substitution rate, along with the many uncertainties 
of the molecular clock 32,33 . Nevertheless, it is worth noting that 
the phylogeographic discontinuity between the N and S clades is 
located close to the northern edge of a mid-peninsular region 
that was repeatedly subjected to marine -flooding during the Plio- 
Pleistocene, and where phylogeographic discontinuities have been 
repeatedly reported 6,28,29,34,35 . Thus, it is plausible that a marine trans- 
gression that occurred in the early Pleistocene might have primed the 
observed divergence. This divergence is conspicuous (sequence 
divergence was 4.0% in the mtDNA data), syntopy was not observed 
between the 2 mtDNA clades, and the admixture of nuclear genomes 
appeared both geographically and quantitatively limited (see 
Figure 3C). Therefore, further study is required to locate the second- 
ary contact zone between both clades and to investigate what micro- 
evolutionary processes occurred following the formation of this zone. 

The transition zone between the NI and Nil lineages in the 
mtDNA data (Figure 2) was not reflected by the transition zone 
between the north-central and the northern population groups in 
the microsatellite data (Figure 3). Although not unusual, it may be 
difficult to reconcile the discordance between nuclear and mitochon- 
drial data to infer the history of populations, given the diversity of 
microevolutionary processes that are potentially involved 36 . 
However, while the northern Apennine is a well known suture zone 
between both intra- and interspecific lineages 27,37-40 , a contact zone in 
the area close to our sites 17-19 would not have any analogues in the 
phylogeographic literature of co-distributed temperate species. In 
addition, no environmental features (present or past) were present 
that could explain the observed pattern. This observation led us to 
speculate that, after the secondary contact between the 2 lineages in 
an area close to the northern Apennines, the NI mitochondrial lin- 
eage incurred massive introgression into the Nil range, resulting in 
the former replacing the latter throughout most of central Italy. On 
the other hand, unfavourable climatic conditions along the central 
Apennines chain during the glacial stages 41,42 might provide a viable 
explanation for the divergence between the 2 lineages. During these 
stages, the NI and Nil lineages might have persisted in temperate 
lowland areas along the eastern Adriatic and the central Tyrrhenian 
coasts, respectively 42 . 

The phylogeographic discontinuity observed between the SI and 
SII mtDNA clades was located in central Calabria. Although a more 
admixed picture emerged at the nuclear microsateilites, central 
Calabria was also the main area of transition between the south- 
central and the southern groups of populations (marked as purple 
and green, respectively, in Figure 3). Interestingly, several phylogeo- 
graphic breaks and secondary contact zones were clustered in this 
area 6,21,27 30 , which was cyclically inundated by marine transgressions 
until the late-Middle/Late Pleistocene 43 . Therefore, it seems plausible 
that these transgressions also primed the divergence between the 2 R. 
dalmatina groups found in south-central and southern Italy. 

Despite the relative geographic contiguity of the alleged past 
ranges of the divergent R. dalmatina lineages found, the lineages 
showed substantially different demographic trends during the Late 
Pleistocene. The star-like structure of the haplotype genealogy for 
lineage NI, its BSP reconstruction, and the estimated values of the 
neutrality test statistics jointly indicated that, after a marked demo- 
graphic bottleneck occurred at the last glacial maximum (LGM; as 
revealed by the estimated TMRCA around 23 000 years ago), NI 



underwent a substantial and abrupt demographic expansion. Since 
NI corresponds to the widespread lineage found by Vences and 
colleagues 26 (2013; data not shown) throughout most of R. dalmati- 
na's range in Europe, and in light of its Middle-Pleistocene diver- 
gence from Nil (0.38 [0.17-0.59] My), we suggest that lineage NI 
re-colonized most of its current range after the LGM from a glacial 
refugium located somewhere between the western Balkans and the 
wide coastal plain that occupied what is now the Adriatic Sea 42,44 . In 
contrast, lineages SI and SII showed the genetic imprints of demo- 
graphic stability, lasting for most of the Late Pleistocene (SI), or 
before (SII). Similar trends of prolonged demographic stability have 
also been hypothesized for intraspecific lineages of other temperate 
species studied to date in southern Italy 6,2 '. Thus, contrary to what is 
frequently observed farther north, the most striking effect of glacial- 
interglacial cycles on the history of populations from south-central 
and southern Italy seems to be driven not by climate forcing on 
range/demographic size variations, but by glacio-eustatic oscillations 
in the sea level and the resulting allopatric fragmentation of popula- 
tions imposed by the intervening seaways. Consequently, this area 
could be viewed as a 'glacial refugium' ('southern refugium' sensu 45 ) 
from the worsening climate during the ice ages and as a set of 'inter- 
glacial refugia' from the rising sea levels during interglacial epochs. 

Among wide-ranging species in Western Europe, geographic pat- 
terns of genetic diversity as observed for R. dalmatina (more diversity 
to the south and less diversity to the north, i.e. the 'southern-richness, 
northern-purity' pattern 16 ) are not unusual 4,7,46 . In most of these 
species (including _R. dalmatina 26 ), the 'southern-richness' compon- 
ent of this pattern, i.e. the existence and location of hotspots, has been 
interpreted as the outcome of the long-term persistence of ancestral 
populations within southern peninsulas 2,3 . In contrast with this view 
for wide-ranging species and in parallel with recent observations in 
species endemic to these peninsulas 5,6,8,47 , our results emphasize the 
key role played by the 'diversity-enhancing' effects of the cycles of 
allopatric divergence within these areas. Since a sampling scheme 
specifically focused on the putative refugial ranges is crucial for dis- 
entangling the above scenarios, we anticipate that as more studies 
adopt such sampling schemes, our view of the diversity of processes 
involved in the formation of hotspots of intraspecific biodiversity will 
continue to change. 

Finally, our results have at least 3 major implications for the con- 
servation of R. dalmatina populations, both within our study area 
and at the whole species' range. First, in order to preserve the species' 
resources for coping with future natural and anthropogenic envir- 
onmental changes, peninsular Italy should receive special attention 
since most of _R. dalmatina 's genetic diversity is apportioned there, 
and in light of the historical role played by this area in the species' 
survival. This appears especially important when considering the 
declining trend recently observed for R. dalmatina populations from 
several portions of its range, including the Italian peninsula 25 . 
Second, since the deep genetic structure of the populations indicates 
the occurrence of at least two ESUs (sensu 23 ) within this area, specific 
conservation efforts should be developed to optimize the amount of 
genetic variation preserved. Finally, a study of the genetic structure of 
the presumed secondary contact zone between the N and S lineages 
in central Italy should be carried out to define the patterns of past and 
present interaction among them, their most appropriate taxonomic 
arrangement, and the best strategy for their conservation. 

Methods 

Sampling and laboratory procedures. We collected tissue samples from 244 
individuals of R. dalmatina from 40 localities, 33 of which were located along the 
Italian peninsula south of the Alps, 3 from Slovenia, and 4 from Croatia. The 
geographical position and number of individuals sampled at each location are given in 
Table 1 and Figure 1. Tissue samples were obtained through toe-clipping, and all 
individuals were released at the place of collection immediately after tissue removal. 
Samples were transported to the laboratory in liquid nitrogen and were stored at 
— 80"C until subsequent analyses. 
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We extracted genomic DNA by using the standard phenol-chloroform protocol 48 . 
Two partial fragments of the mitochondrial DNA (mtDNA) genes encoding CYTB 
and COI were amplified and sequenced. The same PCR protocol was used for both 
mtDNA fragments. Amplifications were performed in a fmal volume of 25 ul con- 
taining MgCl 2 (2.5 raM), reaction buffer (IX; Promega), 4 dNTPs (0.2 mM each), 2 
primers (0.2 uM each), Taq polymerase (2-unit; Go-Taq Promega), and 2 ul of 
template DNA. PCR cycles were as follow: 94 L, C for 5 min followed by 36 cycles of 
94°C for 45 s, annealing temperature (50°C for CYTB and 52°C for COI) for 45 s, 
72'C for 90 s, and a final step at 72"C for 10 min. The primers used for amplification 
and sequencing of the CYTB fragment were 494 49 and MVZ16 50 ; primers for the COI 
fragment were Amp-P3r and Amp-P3f 51 . Sequencing reactions were carried out by 
Macrogen Inc. (www.macrogen.com) by using an ABI PRISM 377 DNA sequencer 
(PE Applied Biosystems) following the ABI PRISM BigDye Terminator Cycle 
Sequencing protocol. All individuals analysed were double sequenced. 

Seven microsatellite loci (F5, G12, Gil, E8, HI, B5, and F3) were amplified using 
previously published protocols 52 . One locus (F3) was subsequently removed from the 
analysis because a large percentage of data were missing (>50%). 

Fragments were sized using the GeneScan 400HD internal size standard, and alleles 
were scored using PEAKSCANNER 1.0 (Applied Biosystems). 

We used the method implemented in the software MICROCHECKER 53 to check 
for the presence of null alleles and large allele dropout; however, there was no evid- 
ence of such problems. 

Mitochondrial DNA data analysis. Sequence electropherograms were checked by 
eye, by using ChromasPro (Technelysium Ltd.), and multiple sequence alignments 
were carried out using ClustalX 54 . Sequences of the CYTB and COI fragments were 
concatenated using Concatenator 1.1. 0 55 , after a partition homogeneity test 56 
conducted in Paup*4.0bl0 did not reject the hypothesis of homogeneity in the 
phylogenetic signal carried by both fragments. The optimal model of sequence 
evolution for both fragments separately and the concatenated dataset were selected 
among 88 possible models using Jmodeltest 0.1.1 57 . The best-fit models selected using 
the Akaike information criterion were TrN + T 58 (gamma distribution shape 
parameter — 0.032) for COI, TrN + I (proportion of invariable sites — 0.073) for 
CYTB, and TrN + I (proportion of invariable sites — 0.824) for the concatenated 
dataset. Therefore, the latter was used in all downstream analyses on the 
unpartitioned dataset. 

Sequences were collapsed into unique haplotypes and their nucleotide variation 
was evaluated with the software DnaSP 5 59 . The same software was used to estimate 
the haplotype (h) and nucleotide (n) diversity 60 among major groups of haplotypes as 
identified by phylogenetic analyses; Mega 5 61 was then used to compute the net 
sequence divergence among haplotype groups. 

The genealogical relationships among the identified haplotypes were investigated 
using 2 approaches. First, we estimated phylogenetic networks using the statistical 
parsimony procedure of Templeton 62 , implemented in the software TCS 1.21 63 using 
the 95% limit for a parsimonious connection. Second, to account for recent concerns 
about the accuracy of network- building methods 64 , we further investigated the 
genealogical relationships among haplotypes by using the approach suggested by 64 . A 
phylogenetic tree topology was estimated using the maximum likelihood algorithm as 
implemented in PhyML 65 . We used the default setting in PhyML for all parameters 
except the substitution model (TrN + I) and the type of tree improvement {SPR & 
NNI 66 ). The estimated tree topology was then converted into a haplotype genealogy 
by using the software Haplotype Viewer 65 . 

Estimates of the divergence time between identified haplogroups were obtained 
using BEAST 1.7.5 67 using the Bayesian Skyline as coalescent tree prior, which 
imposes minimal assumptions on the data 67 . Preliminary analyses were run using a 
relaxed molecular clock model, assuming an uncorrelated lognormal distribution. 
However, since the standard deviation of the uncorrelated lognormal relaxed clock 
(the parameter ucld.stdev) was close to zero, subsequent analyses were run enforcing 
a strict clock model. Since we had neither internal calibration points nor previous 
estimates of the substitution rates for the 2 mtDNA fragments studied in R. dalma- 
Hna, this and subsequent demographic analyses were carried out using a uniform 
prior for the substitution rate, bounded at 0.012 and 0.014 substitutions/site/My, to 
reflect the range of substitution rates previously estimated for CYTB and COI frag- 
ments in other European brown frogs 68,69 . The final analysis was run in triplicate to 
check for consistency among runs, each with a Markov chain Monte Carlo (MCMC) 
run for 20 million generations and sampled every 2000 generations. Further, sat- 
isfactorily high ESS values for all parameters in the analysis were obtained using these 
settings. Convergence, stationarity, and the appropriate number of steps to be dis- 
carded as burn-in were assessed using TRACER 1.6. 

A hierarchical analysis of molecular variance (AMOVA 70 ) was carried out by using 
ARLEQUIN 3.1 71 , in order to estimate the amount of variation attributable to dif- 
ferences among population groups, among populations within groups, and within 
populations. This analysis was run using the TrN model, with population groups 
defined according to the geographical distribution of the main haplogroups. 

The historical demography of the main haplogroups was inferred using the 
Bayesian skyline plot method (BSP) implemented in BEAST 72,73 . Five preliminary 
replicates of the BSP analysis were carried out to fme-tune the parameters of the 
MCMC using the auto -optimization option. Then, we ran the final analyses in tri- 
plicates for each haplogroup to check for consistency. The number of groups (m) was 
set to 10 (but setting different values in preliminary analyses did not affect the 
resulting demographic trends); the MCMCs were run for 30 million generations and 
sampled every 1000 generations, with the first 10% discarded as burn-in. BSP analyses 



were run using the TrN + I, a strict clock model, and a uniform prior for the 
substitution rate, bounded at 0.012 and 0.014 substitutions/site/My. The results were 
analysed and summarised using TRACER 1.6. Finally, the occurrence of past 
demographic changes was evaluated by computing the neutrality test statistics F s 74 , 
R 2 75 and D 76 , by using DNaSP. The significance of these statistics was assessed using 
1000 coalescent simulations. 

Microsatellites data analysis. We tested for deviations from the expected Hardy- 
Weinberg and linkage equilibria using FSTAT 2.9.3. The software GENETIX was 
used to obtain estimates of genetic diversity based on the observed and the expected 
heterozygosity and the number of alleles per locus. 

The extent of population genetic structure across the study area was investigated by 
the Bayesian clustering algorithm implemented in TESS 2.3 77 . This method is 
advantageous in that it takes the spatial distribution of sampled individuals into 
account, and has been shown to outperform similar methods, particularly in cases of 
shallow population structures 77 . The TESS analysis was carried out by modelling 
admixture using the conditional autoregressive (CAR) model. A preliminary analysis 
with 20 000 runs (of which the first 5 000 were discarded as burn-in) and 10 replicates 
for each K value (i.e. the number of clusters) between 2 and 20, was carried to assess 
model performance, the need to fme-tune the model parameters, and to reduce the 
range of plausible K values. The fmal analysis contained 100 replicates with K = 2 to 7, 
each with 50 000 runs (of which the first 20 000 were discarded as burn-in). The 
spatial interaction parameter was set to 0.6 (i.e. the default value) and the option to 
update this parameter was activated. The model that best fit the data was selected 
using the deviance information criterion (DIC). DIC values were averaged over the 
100 replicates for each K value, and the most likely K value was selected as the one at 
which the average DIC reached a plateau. The estimated admixture proportions of the 
10 runs with the lowest DIC values were averaged using CLUMPP 1.1. 2 78 . Finally, to 
obtain an easily interpretable representation of the estimated admixture proportions, 
we generated posterior predictive maps by the spatial interpolation procedure 
available in TESS, and with the maps being visualised by TESS Ad-Mixer 79 . Finally, we 
used the posterior estimate of the allele frequencies at each locus in each cluster to 
compute the mean expected heterozygosity of each cluster. 
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